% Example on Multiplication Quantization effects on the 
%  First-Order IIR filter:
%   y(n) = Q[-a*y(n-1)] + x(n); 
%   x(n) scaled to prevent overflow

% Example Parameters
B = 10;               % # of fractional bits
BB = 2^B;             % useful factor in quantization
N = 10000;            % # of samples
xn = (2*rand(1,N)-1); % Input sequence - Uniform Distribution
a = 0.9;              % Filter parameter
Xm = 1-abs(a);        % Scaling factor


% Quantize the input and the filter coefficients
  xn = round(Xm*xn*BB)/BB; % Scaled Input quantized to B bits
   a = round(a*BB)/BB;     % a quantized to B bits

% Filter output without multiplication quantization
  yn = filter(1,[1,-a],xn);  % output using filter routine

% Filter output with multiplication quantization
   yq = zeros(1,N);                        % Initialize quantized output array
yq(1) = round(b(1)*xt(1)*DD)/DD;           % Calculation of the first sample yq(1)
  B0X = round(b(1)*xt(2)*DD)/DD;           % Quantization of b0*x(2)
  B1X = round(b(2)*xt(1)*DD)/DD;           % Quantization of b1*x(1)
  A1Y = round(a(1)*yq(1)*DD)/DD;           % Quantization of a1*y(1)
yq(2) = B0X+B1X+A1Y;                       % Calculation of the second sample yq(2)
for I = 3:N;
    B0X = round(b(1)*xt(I)*DD)/DD;         % Quantization of b0*x(n)
    B1X = round(b(2)*xt(I-1)*DD)/DD;       % Quantization of b1*x(n-1)
    B2X = round(b(3)*xt(I-2)*DD)/DD;       % Quantization of b2*x(n-2)
    A1Y = round(a(1)*yq(I-1)*DD)/DD;       % Quantization of a1*y(n-1)
    A2Y = round(a(2)*yq(I-2)*DD)/DD;       % Quantization of a2*y(n-2)
  yq(I) = B0X+B1X+B2X+A1Y+A2Y;             % Calculation of the I-th sample yq(I)
end

% Output Error Analysis
   en = yt-yq; clear yt yq;                % Outout error sequence
eemax = max(en); eemin = min(en);          % Maximun and minimum of the error
enmax = max(abs([eemax,eemin]));           % Absolute maximum range of the error
enavg = mean(en); enstd = std(en);         % Mean and std dev of the error
   en = round(en*(10^bM)/(2*enmax)+0.5);   % Normalized en (interger between -M & M)
   en = sort([en,-M:1:(M+1)]);             % 
    H = diff(find(diff(en)))-1; clear en;  % Error histogram
    H = H/N;                               % Normalized histogram
 Hmax = max(H); Hmin = min(H);             % Max and Min of the normalized histogram

% Plot variables
  bM = 2; DbM = 2^bM;                      % bin parameter
   M = round(DbM/2);                       % Half number of bins
bins = [-M+0.5:1:M-0.5];                   % Bin values from -M to M
   Q = bins/DbM;                           % Normalized bins
 YTN = 10^(-bM);                           % Ytick marks interval
 YLM = 4.0*YTN;                            % Yaxis limit

% Plot of the error histogram
Hf_1 = figure('Units','normalized','position',[0.05,0.05,0.9,0.9],'color',[0,0,0]);
set(Hf_1,'NumberTitle','off','Name','Experiment 4-3 : Direct Form-1');

bar(Q,H); axis([-0.5,0.5,0,YLM]);
title(TITLE,'FontSize',16,'FontName','Times', 'Fontweight','Bold');
set(gca,'TickDir','in','FontWeight','Bold'); 
xlabel('Normalized error'); ylabel('Distribution of error');
set(gca,'XTickMode','manual','XTick',[-0.5:0.1:0.5]);
set(gca,'YTickmode','manual','YTick',[0:YTN:YLM]);
LABEL1 = ['SAMPLE SIZE N = ', num2str(N)];
LABEL2 = ['ROUNDED T0 D  = ', num2str(D), ' DECIMALS'];
LABEL3 = ['MEAN * 10E+D  = ', num2str(enavg*DD)];
LABEL4 = ['SIGMA * 10E+D = ', num2str(enstd*DD)];
LABEL5 = ['MIN PROB BAR HEIGHT = ', num2str(Hmin)];
LABEL6 = ['MAX PROB BAR HEIGHT = ', num2str(Hmax)];
LABEL7 = ['MIN ERROR * 10E+D   = ', num2str(eemin*DD)];
LABEL8 = ['MAX ERROR * 10E+D   = ', num2str(eemax*DD)];
text(-0.47,3.80*YTN,LABEL1,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(-0.47,3.65*YTN,LABEL2,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(-0.47,3.50*YTN,LABEL3,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(-0.47,3.35*YTN,LABEL4,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.80*YTN,LABEL5,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.65*YTN,LABEL6,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.50*YTN,LABEL7,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.35*YTN,LABEL8,'Fontname','Courier','FontWeight','Bold','Fontsize',10);